#ifndef __CURVENDIMENSION_HPP__

#define __CURVENDIMENSION_HPP__


#include <vector>
#include <algorithm>
#include "util.h"
#include <QDebug>


template < typename real , unsigned int dimen , class Point , class Vector >
class CurvePoligonal {
public:
    CurvePoligonal( void ) { pIsClosed = false ; }

    CurvePoligonal( const CurvePoligonal<real,dimen,Point,Vector> &curve ) { this->copy(curve) ; }

    CurvePoligonal( const std::vector<Point> &points ){
        pPoints = points;
        if ( pPoints.size() != 0 && points[ pPoints.size()-1 ] == pPoints[ 0 ] )
        {
            pIsClosed = true ;
            pPoints.pop_back();
        }
        else pIsClosed = false ;
    }

    virtual void copy( const CurvePoligonal<real,dimen,Point,Vector> &curve ) {
        pPoints = curve.pPoints;
        pIsClosed = curve.pIsClosed;
    }

    virtual CurvePoligonal &operator=( const CurvePoligonal<real,dimen,Point,Vector> &curve ) {
        this->copy(curve);
        return (*this);
    }

    virtual Point& operator[](int i) {
        unsigned idx = i;
        if(i < 0)
            idx = i + pPoints.size();
        if(i >= (int)pPoints.size())
            idx = i - pPoints.size();

        return pPoints[idx] ;
    }
    virtual Point  operator[](int i) const {
        unsigned idx = i;
        if(i < 0)
            idx = i + pPoints.size();
        if(i >= (int)pPoints.size())
            idx = i - pPoints.size();
        return pPoints[idx] ;
    }

    // Number of points of the curve
    virtual unsigned int size( void ) const { return pPoints.size() ; }

    virtual void setPoints( const std::vector<Point> &points ) { pPoints = points ; }
    virtual std::vector<Point> getPoints( void ) const { return pPoints ; }

    virtual void clear( void ) { pPoints.clear() ; pIsClosed = false; }

    virtual Point   at( unsigned int i ) const {
        AssertBound(i<size()) ;
        return pPoints[i] ;
    }
    virtual Point   atBegin( void ) const { AssertBound(size()>0) ; return pPoints[0] ; }
    virtual Point   atEnd( void ) const {  AssertBound(size()>0) ; return pPoints[ size()-1 ] ; }

    virtual Point   eval( real u ) const { Point p ; eval( u , p ) ; return p ; }
    virtual unsigned int  eval( real u , Point& p ) const { real alpha ; return eval( u , p , alpha ); }
    // return the index of the point before u be archived, p is the point on the curve in u and
    // alpha is the relative distance of the point at(index) to the point p.
    // It means if alpha = 0 then p = at(index) and if alpha = 1 then p = at(index+1)
    virtual unsigned int  eval( real u , Point& p , real& alpha ) const {
        if( size() == 0 ) return 0;
        real l = length();
        if( u >= 1.0 )
        {
            p = at(size()-1);
            return size()-1;
        }
        if( u <= 0.0 || l == 0 )
        {
            p = at(0);
            alpha = 0.0;
            return 0;
        }
        real il = 1.0/l;
        real aux1 = 0 ;
        real aux2 = 0 ;
        real lenTmp = 0 ;
        unsigned int index = 0 ;
        Vector vtmp;
        while( aux1 <= u )
        {
            lenTmp = aux1;
            index++;
            vtmp = pPoints[index]-pPoints[index-1];
            aux2 = vtmp.norm()*il;
            aux1 += aux2;
        }
        alpha = (u - lenTmp)/aux2;
        index--;
        p = at(index) + alpha*vtmp;
        return index;
    }

    virtual void setPoint( unsigned int index , const Point& p ){ AssertBound(index<size()) ; pPoints[index] = p ; }
    virtual void add( const Point& p ) { pPoints.push_back(p) ;}

    virtual void insert( unsigned int index , const Point& p ) { if( index < size()) pPoints.insert( pPoints.begin()+index , p ); else add(p) ; }
    virtual void remove(int index) { unsigned int idx = index%size(); pPoints.erase(pPoints.begin()+idx); }
    virtual int remove(std::vector<int> list)
    {
        for(unsigned int i = 0; i < list.size(); ++i )
        {
            if(list[i] < 0)
                list[i] = list[i] + pPoints.size();
            if(list[i] >= (int)pPoints.size())
                list[i] = list[i] - pPoints.size();
        }

        //qDebug() << trans(list);

        std::sort(list.begin(), list.end());
        for(int i = 0; i < list.size()/2; ++i)
        {
            int t = list[i];
            list[i] = list[list.size() -i -1];
            list[list.size() -i -1] = t;
        }

        //qDebug() << trans(list);

        for(unsigned int i = 0; i < list.size(); ++i )
            remove(list[i]);

        return list[list.size()-1];
    }

    virtual void pop_back(void) { pPoints.pop_back(); }
    virtual void push_back( const Point& p ) { pPoints.push_back(p) ; }

    virtual void close( bool b = true ) { pIsClosed = b ; }
    virtual void open( void ) { pIsClosed = false ; }
    virtual bool isClosed() const { return pIsClosed ; }

    virtual void chaikinFilter( unsigned int numberOfSimplifications = 1 ){
        if (size() == 0) return;
        if( isClosed() )
        {
            for( unsigned int n = 0 ; n < numberOfSimplifications ; ++n )
            {
                std::vector<Point> tmp;
                // make sure we do not loose any points!!
                if ( size() <= 6 ) return;

                // step over every 2 points
                for( unsigned int i = 1 ; i < ( size() - 2 ) ; i += 2 )
                {
                    // get original points
                    const Vector p0 = pPoints[i-1].toVector();
                    const Vector p1 = pPoints[i  ].toVector();
                    const Vector p2 = pPoints[i+1].toVector();
                    const Vector p3 = pPoints[i+2].toVector();

                    // calculate the original point
                    Vector Q = -0.25f*p0 + 0.75f*p1 + 0.75f*p2 - 0.25f*p3;

                    // add to new curve
                    tmp.push_back( Point(Q) );
                }
                if( size()%2 != 0 )
                {
                    unsigned int i = size() - 1 ;
                    // get original points
                    const Vector p0 = pPoints[i-1].toVector();
                    const Vector p1 = pPoints[i  ].toVector();
                    const Vector p2 = pPoints[0].toVector();
                    const Vector p3 = pPoints[1].toVector();
                    // calculate the original point
                    Vector Q = -0.25f*p0 + 0.75f*p1 + 0.75f*p2 - 0.25f*p3;

                    // add to new curve
                    tmp.push_back( Point(Q) );
                }
                else
                {
                    unsigned int i = size() - 2 ;
                    // get original points
                    const Vector p0 = pPoints[i-1].toVector();
                    const Vector p1 = pPoints[i  ].toVector();
                    const Vector p2 = pPoints[i+1].toVector();
                    const Vector p3 = pPoints[0].toVector();
                    // calculate the original point
                    Vector Q = -0.25f*p0 + 0.75f*p1 + 0.75f*p2 - 0.25f*p3;

                    // add to new curve
                    tmp.push_back( Point(Q) );

                }
                // copy over pPoints
                pPoints = tmp ;
            }
        }
        else
        {
            for( unsigned int n = 0 ; n < numberOfSimplifications ; ++n )
            {
                std::vector<Point>  tmp;
                // make sure we do not loose any points!!
                if (size() <= 8) return;

                // keep the first point
                tmp.push_back( pPoints[0] );
                Point p2 = -0.5f*pPoints[0].toVector() + pPoints[1].toVector() +
                        0.75f*pPoints[2].toVector() - 0.25f*pPoints[3].toVector();
                tmp.push_back( Point(p2) );

                // step over every 2 points
                for( unsigned int i = 2 ; i < ( size() - 5 ) ; i += 2 )
                {
                    // get original points
                    const Vector p0 = pPoints[i    ].toVector();
                    const Vector p1 = pPoints[i + 1].toVector();
                    const Vector p2 = pPoints[i + 2].toVector();
                    const Vector p3 = pPoints[i + 3].toVector();
                    // calculate the original point
                    Vector Q = -0.25f*p0 + 0.75f*p1 + 0.75f*p2 - 0.25f*p3 ;
                    // add to new curve
                    tmp.push_back( Point(Q) );
                }
                unsigned int lastIndex = size() - 1 ;
                Vector pL = -0.25f*pPoints[lastIndex-3].toVector() + 0.75f*pPoints[lastIndex-2].toVector() +
                        pPoints[lastIndex-1].toVector() - 0.50f*pPoints[lastIndex].toVector();
                tmp.push_back( Point(pL) );
                tmp.push_back(pPoints[ lastIndex ]);
                // copy over pPoints
                pPoints = tmp ;
            }
        }
    }

    virtual void chaikinSubDivide( unsigned int numberOfsimplifications = 1 ){
        if (size() == 0) return;
        if( isClosed() )
        {
            for( unsigned int n = 0 ; n < numberOfsimplifications ; ++n )
            {
                std::vector<Point>  tmp;
                // step over every 2 points
                for( unsigned int i = 1 ; i < size() ; ++i )
                {
                    // get original points
                    const Vector p0 = pPoints[i-1].toVector();
                    const Vector p1 = pPoints[i  ].toVector();

                    // calculate the original point
                    Vector Q = 0.75f*p0 + 0.25f*p1 ;

                    tmp.push_back( Point(Q) );

                    Vector R ;
                    R = 0.25f*p0 + 0.75f*p1 ;
                    tmp.push_back( Point(R) );
                }
                // get original points
                const Vector p0 = atEnd().toVector();
                const Vector p1 = atBegin().toVector();

                // calculate the original point
                Vector Q = 0.75f*p0 + 0.25f*p1 ;

                tmp.push_back( Point(Q) );

                Vector R ;
                R = 0.25f*p0 + 0.75f*p1 ;
                tmp.push_back( Point(R) );
                // copy over pPoints
                pPoints = tmp ;
            }
        }
        else
        {
            for( unsigned int n = 0 ; n < numberOfsimplifications ; ++n )
            {
                std::vector<Point>  tmp;
                // keep the first point
                tmp.push_back( pPoints[0] );

                // step over every 2 points
                for( unsigned int i = 1 ; i < size() ; ++i )
                {
                    // get original points
                    const Vector p0 = pPoints[i-1].toVector();
                    const Vector p1 = pPoints[i  ].toVector();

                    // calculate the original point
                    Vector Q = 0.75f*p0 + 0.25f*p1 ;
                    tmp.push_back( Point(Q) );
                    Vector R ;
                    R = 0.25f*p0 + 0.75f*p1 ;
                    tmp.push_back( Point(R) );
                }
                tmp.push_back(pPoints[ size()-1 ] );
                // copy over pPoints
                pPoints = tmp ;
            }
        }
    }

    virtual void superSample( real step ){
        if (size() == 0) return;

        CurvePoligonal tmpLine;
        tmpLine.add( pPoints[0] );
        for( unsigned int i = 0 ; i <  size() - 1 ; ++i )
        {
            Vector tmpV = pPoints[i+1] -pPoints[i] ;
            real tmpNorm = tmpV.norm() ;
            if( tmpNorm > step )
            {
                Vector base = pPoints[i].toVector() ;
                real invNorm = 1/tmpNorm ;
                Vector unitV = invNorm*tmpV ;
                unsigned int j = 1 ;
                while( (j*step) <= ( tmpNorm - step ) )
                {
                    Vector newP = base + (j*step)*unitV ;
                    tmpLine.add( Point( newP ) );
                    ++j;
                }
            }
            tmpV = pPoints[i+1] -tmpLine.atEnd() ;
            tmpNorm = tmpV.norm() ;
            if( tmpNorm > 0.8*step )
            {
                tmpLine.add( pPoints[i+1] );
            }
        }
        if(isClosed())
        {
            Vector tmpV = atBegin() - atEnd() ;
            real tmpNorm = norm( tmpV ) ;
            //        real step = 0.5f;
            if( tmpNorm > step )
            {
                Vector base = atEnd().toVector() ;
                real invNorm = 1/tmpNorm ;
                Vector unitV = invNorm*tmpV ;
                unsigned int j = 1 ;
                while( (j*step) <= ( tmpNorm - step ) )
                {
                    Vector newP = base + (j*step)*unitV ;
                    tmpLine.add( Point( newP ) );
                    ++j;
                }
            }
        }
        tmpLine.close( isClosed() );
        *this = tmpLine;
    }

    // apply the superSample(step) and chaikinFilter( numberOfInteractions )
    // usefull to cleanup hand-draw input curves
    virtual void lineFilter( real step = 0.5 , unsigned int numberOfInteractions = 5 ){
        if (size() == 0) return;
        this->superSample( step );
        this->chaikinFilter( numberOfInteractions );
    }

    virtual void meanFilter( void ){
        if( size() > 1 )
        {
            Vector v;
            std::vector<Point> tmp;
            if(isClosed())
            {
                v = ( atEnd().toVector()+3*at(0).toVector()+at(1).toVector() )/5.0;
                tmp.push_back( Point(v[0],v[1],v[2]) );
            }
            else tmp.push_back( atBegin() );

            for( unsigned int i = 1 ; i < size()-1 ; ++i )
            {
                v = ( at(i-1).toVector()+3*at(i).toVector()+at(i+1).toVector() )/5.0;
                tmp.push_back( Point(v[0],v[1],v[2]) );
            }
            if (isClosed())
            {
                v = ( at(size()-2).toVector()+3*atEnd().toVector()+at(0).toVector() )/5.0;
                tmp.push_back( Point(v[0],v[1],v[2]) );
            }
            else tmp.push_back( atEnd() );
            pPoints = tmp;
        }
    }

    virtual void reverse( void ){
        std::vector<Point> tmp;
        for( unsigned int i = size() ; i > 0 ; --i )
        {
            tmp.push_back( pPoints[i-1] );
        }
        pPoints = tmp;
    }

    // Join two curves but it test if the first and last points
    // are not repeted, i.e., they are nearer than err.
    virtual void join( const CurvePoligonal& l , real err = EPS ){
        if(l.size() == 0 ) return;
        real dist = INF ;
        if( this->size() > 0 ) dist = ( atEnd() - l.atBegin()).norm();
        unsigned int start = ( dist < err ) ? 1 : 0 ;
        for( unsigned int i = start ; i < l.size()-1 ; ++i ) this->add(l.at(i));
        dist = (l.atEnd() - atBegin()).norm() ;
        if( dist < err ) close( true );
        else add(l.atEnd());
    }

    virtual Vector tangentEval( real u ) const {
        Vector tan;
        if( size() > 1 )
        {
            Point p;
            real alpha;
            unsigned int index = eval( u , p , alpha );
            if( alpha > EPS )
            {
                if( index != size()-1 ) tan = at(index+1) - at(index);
                else if ( isClosed() ) tan = at(0) - at(index);
                else tan = at(index-1) - at(index);
            }
            else tan = tangent(index);
        }
        return tan;
    }

    virtual Vector tangent( unsigned int index ) const {
        Vector tan;
        if( size() > 1 )
        {
            if( index < size()-1 && index > 0 ) tan = ( at(index+1) - at(index-1) )*.5 ;
            else if( isClosed() ) tan = ( at( (index+1)%size() ) - at( (int(index) - 1)%size() ) ) *.5;
            else if ( index > 0 ) tan = at(index) - at(index-1);
            else tan = at(1) - at(0);
        }
        return tan;
    }

    virtual Vector tangentEvalContinuos( real u ) const {
        Vector tan;
        if( size() > 1 )
        {
            Point p;
            real alpha;
            unsigned int index = eval( u , p , alpha );
            if( index != size()-1 ) tan = alpha * tangent(index) + (1.0-alpha)* tangent(index+1);
            else  tan = tangent(index);
        }
        return tan;
    }

    // return the bounding box os the curve
    virtual void minMax( Point& min , Point& max ) const {
        max = Point( -INF ) , min = Point( INF ) ;

        for( unsigned int i = 0 ; i < size() ; i++ )
        {
            for( unsigned int j = 0 ; j < dimen ; j++ )
            {
                max[j] = ( at(i)[j] > max[j] ) ? at(i)[j] : max[j] ;
                min[j] = ( at(i)[j] < min[j] ) ? at(i)[j] : min[j] ;
            }
        }

    }

    virtual real length( void ) const {
        real len = 0;
        Vector vtmp;
        for( unsigned int i = 1 ; i < size() ; i++ )
        {
            vtmp = pPoints[i]-pPoints[i-1];
            len += vtmp.norm();
        }
        if(isClosed()) vtmp = atBegin() - atEnd() , len += vtmp.norm() ;
        return len;
    }
    virtual real length( unsigned int k ) const {
        real len = 0;

        k = ( k < size() )? k+1 : size();

        for( unsigned int i = 1 ; i < k ; i++ )
        {
            Vector vtmp = pPoints[i]-pPoints[i-1];
            len += vtmp.norm();
        }
        return len;
    }

    virtual real getParameter( unsigned int i ) const { return length(i)/length() ; }

    virtual void translate ( const Vector& v ) { for( unsigned int i = 0 ; i < size() ; ++i ) pPoints[i] += v; }

    virtual real distanceTo( const Point& p ) const { unsigned int index ; real u ; Point pr ; return projectPoint( p , index , u , pr ) ; }
    // Project p over this curve and return: the distance between p and pr, its projection, in addition
    // index is the index of the vertex before pr and parameterU is its parameter in the curve
    virtual real projectPoint( const Point& p , unsigned int& index , real& parameterU , Point& pr ) const {
        real dist = INF, uTmp1 = 0.0, uTmp2 = 0.0 ;
        parameterU = 0.0;
        if( size() > 0 )
        {
            index = 0 ;
            Point prTmp = at(0);
            pr = prTmp;
            dist = ( atBegin() - p ).norm() ;
            real distTMP = INF , invL = 0.0 ;
            if( length()>EPS) invL = 1.0/length();
            for( unsigned int i = 1 ; i < size() ; i++ )
            {
                uTmp1 = uTmp2;
                unsigned int iTmp = i;
                Vector u = at(i) - at(i-1);
                Vector v = p - at(i-1);
                Vector w = p - at(i);
                real normQ = u.norm2();
                real alpha = 0.0;
                real distU = u.norm();
                real distV = v.norm();
                real distW = w.norm();

                if( normQ > EPS ) alpha = (u*v)/normQ ;

                if( alpha > 0.0 && alpha < 1.0 )
                {
                    prTmp = at(i-1) + alpha*u;
                    distTMP = (prTmp-p).norm();
                    uTmp1 += distU*invL*alpha;
                }
                else if ( distV < distW ) distTMP =  distV, prTmp = at(i-1) ;
                else distTMP = distW , ++iTmp , uTmp1 += distU*invL , prTmp = at(i) ;

                if( distTMP < dist ) dist = distTMP , index = iTmp - 1 , parameterU = uTmp1 , pr = prTmp;
                uTmp2 += distU*invL;
            }
            if(isClosed())
            {
                uTmp1 = uTmp2;
                unsigned int iTmp =  size();
                Vector u = atBegin() - atEnd();
                Vector v = p - atEnd();
                Vector w = p - atBegin();
                real normQ = u.norm2();
                real alpha = 0.0;
                real distU = u.norm();
                real distV = v.norm();
                real distW = w.norm();

                if( normQ > EPS ) alpha = (u*v)/normQ ;

                if( alpha > 0.0 && alpha < 1.0 )
                {
                    prTmp = atEnd() + alpha*u;
                    distTMP = (prTmp-p).norm();
                    uTmp1 += distU*invL*alpha;
                }
                else if ( distV < distW ) distTMP =  distV, prTmp = atEnd() ;
                else distTMP = distW , iTmp = 0 , uTmp1 += distU*invL , prTmp = atBegin() ;

                if( distTMP < dist ) dist = distTMP , index = iTmp - 1 , parameterU = uTmp1 , pr = prTmp;
            }
        }
        return dist;
    }

    virtual void scale(const Vector & v) {
        for( unsigned int i = 0 ; i < size() ; i++ )
            for( unsigned int j = 0 ; j < dimen ; j++ ) pPoints[i][j]*=v[j];
    }

    virtual Point centroid( void ) const {
        Vector barTmp;

        for( unsigned int i = 0 ; i < size() ; i++ ) barTmp += at(i).toVector();

        barTmp /= size() ;

        return Point( barTmp[0] , barTmp[1] , barTmp[2] );
 }

    // split this curve in the index this curve will be the part before index and returns the curve after,
    // the point at(index) is in both
    virtual CurvePoligonal<real,dimen,Point,Vector> split( unsigned int index ) {

        CurvePoligonal<real,dimen,Point,Vector> newCurve;

        if(index<size())
        {
            std::vector<Point> tmp ;
            tmp = std::vector<Point>( pPoints.begin()+index , pPoints.end() ) ;
            newCurve.setPoints(tmp) ;
            tmp = std::vector<Point>( pPoints.begin(),pPoints.begin()+index+1 ) ;
            setPoints(tmp);
            if(isClosed()) newCurve.add( atBegin() );
            open();
            newCurve.open();
        }

        return newCurve;
    }
    // create 2 or 3 new curves, if this curve is closed: 2 (c1 and c2); if not: 3 (c1, c2 and c3)
    // the new curves are this curve split in the index1 and index2
    virtual bool split( unsigned int index0 , unsigned int index1 , CurvePoligonal<real,dimen,Point,Vector>& c1 ,
                        CurvePoligonal<real,dimen,Point,Vector>& c2 , CurvePoligonal<real,dimen,Point,Vector>& c3  ) const {
        bool b = index1<size() && index0<size() ;
        if( b )
        {
            std::vector<Point> tmp ;
            if(isClosed())
            {
                tmp = std::vector<Point>( pPoints.begin() , pPoints.begin()+index0+1 ) ;
                c1.setPoints(tmp) ;
                tmp = std::vector<Point>( pPoints.begin()+index0 , pPoints.begin()+index1+1 ) ;
                c2.setPoints(tmp) ;
                tmp = std::vector<Point>( pPoints.begin()+index1 , pPoints.end() ) ;
                c3.setPoints(tmp) ;
                c3.join(c1);
                c1 = c3;
                c3.clear();
            } else
            {
                tmp = std::vector<Point>( pPoints.begin() , pPoints.begin()+index0+1 ) ;
                c1.setPoints(tmp) ;
                tmp = std::vector<Point>( pPoints.begin()+index0 , pPoints.begin()+index1+1 ) ;
                c2.setPoints(tmp) ;
                tmp = std::vector<Point>( pPoints.begin()+index1 , pPoints.end() ) ;
                c3.setPoints(tmp) ;
            }
        }
        return b;
    }
    // create 2 new curves (c1 and c2) the new curves are this curve split in the index
    virtual bool split( unsigned int index , CurvePoligonal<real,dimen,Point,Vector>& c1 ,
                        CurvePoligonal<real,dimen,Point,Vector>& c2 ) const {
        c1 = (*this);
        c2 = c1.split( index );
        return c2.size() > 0;
    }
    // create 2 new curves (c1 and c2) the new curves are this curve split in the projection of p, if addPoint is true,
    // p is add to c1 and c2 instead of pr
    virtual bool split( const Point& p , CurvePoligonal<real,dimen,Point,Vector>& c1 ,
                        CurvePoligonal<real,dimen,Point,Vector>& c2 , bool addPoint = false , real err = EPS ) const {
        real u;
        unsigned int index;
        Point pr;

        c1 = (*this);

        c1.projectPoint( p , index , u , pr );
        index += (index>0) ? 0 : 1 ;
        if(addPoint)
        {
            if( (c1.at(index)-p).norm() > err && (c1.at(index-1)-p).norm() > err )
                ++index , c1.insert(index,p);
        }
        else if( (c1.at(index)-pr).norm() > err && (c1.at(index-1)-pr).norm() > err )
        {
            ++index ;
            std::cerr << " index = " << index << std::endl;
            c1.insert( index , pr );
        }
        c2 = c1.split(index) ;
        return c2.size() > 0;
    }

    // return a new curve that is the result of use of this curve to over-sketch curve, the difference between
    // the returned curve and curve is saved in rest, this curve could have the orientation changed.
    // The parameter "err" controls if two points are the same and aproximationFactor is how many times
    // the start and or points should be far of the curve to be interpreted as extension not as over-sketch
    virtual CurvePoligonal<real,dimen,Point,Vector> overSketch( const CurvePoligonal<real,dimen,Point,Vector>& curve ,
                                                                CurvePoligonal<real,dimen,Point,Vector>& rest ,
                                                                real err = EPS , real aproximationFactor = 64.0 ) {
        Point p0 = atBegin();
        Point p1 = atEnd();

        Point pr0 , pr1 , prTmp;
        real u0, u1, uTmp, distTmp;
        unsigned int index0 , index1, indexTmp;

        real dist0 = curve.projectPoint( p0 , index0 , u0 , pr0 ) ;
        real dist1 = curve.projectPoint( p1 , index1 , u1 , pr1 ) ;

        if( u0 > u1 )
        {
            prTmp = pr0;
            uTmp = u0;
            indexTmp = index0;
            distTmp = dist0;
            pr0 = pr1;
            u0 = u1;
            index0 = index1;
            dist0 = dist1;
            pr1 = prTmp;
            u1 = uTmp;
            index1 = indexTmp;
            dist1 = distTmp;
            this->reverse();
        }

        CurvePoligonal<real,dimen,Point,Vector> A , B , C , curveTmp0 , curveTmp1  ;

        if( dist0 > aproximationFactor*err && dist1 < aproximationFactor*err )
        {
            curve.split( pr1 , rest , curveTmp0 , false , err );
            this->join( curveTmp0 , err );
            curveTmp0 = (*this);
        }
        else if( dist0 < aproximationFactor*err && dist1 > aproximationFactor*err   )
        {
            curve.split( pr0 , curveTmp0 , rest , false , err );
            curveTmp0.join( (*this) , err ) ;
        }
        else
        {
            curve.split( pr0 , A , curveTmp0 , false , err );
            curveTmp0.split( pr1 , B , C , false , err );

            Point centOr = curve.centroid();

            curveTmp0 = A;
            curveTmp0.join((*this),err);
            curveTmp0.join(C,err);
            curveTmp0.close(isClosed());
            Point cent0 = curveTmp0.centroid();

            curveTmp1 = B;
            curveTmp1.reverse();
            curveTmp1.join((*this),err);
            curveTmp1.close();
            Point cent1 = curveTmp1.centroid();

            if( (centOr-cent0).norm2() < (centOr-cent1).norm2() ) rest = B ;
            else
            {
                curveTmp0 = curveTmp1;
                rest = C;
                rest.join(A,err);
            }
        }

        return curveTmp0 ;
    }

private:
    std::vector<Point> pPoints;
    bool pIsClosed;

};


#endif // __CURVENDIMENSION_HPP__

